Load packages

library(tidyverse)
library(emmeans)
library(glmmTMB)
library(kableExtra)
library(lme4)
library(lmerTest)

Load dataset

data <- read_csv("SB12.csv")

Glimpse dataset

glimpse(data)

Rows: 160
Columns: 19
$ location                 <chr> "Lancaster", "Lancaster", "Lancaster", "Lanc…
$ year                     <dbl> 2019, 2019, 2019, 2019, 2019, 2019, 2019, 20…
$ herb                     <chr> "PRE only", "PRE only", "PRE only", "PRE onl…
$ other                    <chr> "PRE only", "PRE only", "PRE only", "PRE onl…
$ trade                    <chr> "PRE only", "PRE only", "PRE only", "PRE onl…
$ plot                     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ trt                      <dbl> 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4,…
$ rep                      <dbl> 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3,…
$ waterhempcontrol_14      <dbl> 85, 85, 95, 35, 99, 99, 99, 99, 91, 95, 85, …
$ grasscontrol_14          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ waterhempcontrol_28      <dbl> 88, 25, 95, 10, 95, 99, 99, 99, 98, 88, 75, …
$ grasscontrol_28          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ waterhempcontrol_harvest <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ counts_m2                <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ biomass_gm2              <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ yield_bu                 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ phyto_14                 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ dicambadrift_14          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ canopyclosure_14         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …

Skim dataset

skimr::skim(data)

Data summary

Name	data
Number of rows	160
Number of columns	19
_______________________
Column type frequency:
character	4
numeric	15
________________________
Group variables	None

Variable type: character

skim_variable	complete_rate	min	max	n_unique
location	1	8	9	2
herb	1	8	46	10
other	1	8	46	10
trade	1	5	14	10

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
year	0	1.00	2019.50	0.50	2019.00	2019.00	2019.50	2020.00	2020.00	▇▁▁▁▇
plot	80	0.50	255.50	112.55	101.00	178.25	255.50	332.75	410.00	▇▇▁▇▇
trt	0	1.00	5.50	2.88	1.00	3.00	5.50	8.00	10.00	▇▇▇▇▇
rep	0	1.00	2.50	1.12	1.00	1.75	2.50	3.25	4.00	▇▇▁▇▇
waterhempcontrol_14	0	1.00	92.71	12.59	35.00	92.00	99.00	99.00	100.00	▁▁▁▁▇
grasscontrol_14	124	0.22	87.36	10.73	65.00	79.00	89.50	98.00	100.00	▂▃▃▃▇
waterhempcontrol_28	0	1.00	85.47	19.62	0.00	80.00	94.00	99.00	100.00	▁▁▁▂▇
grasscontrol_28	120	0.25	90.65	25.88	0.00	100.00	100.00	100.00	100.00	▁▁▁▁▇
waterhempcontrol_harvest	80	0.50	78.47	28.66	0.00	75.00	88.50	95.00	99.00	▁▁▁▂▇
counts_m2	80	0.50	3.05	3.55	0.00	1.00	2.00	4.00	18.00	▇▂▁▁▁
biomass_gm2	80	0.50	17.82	51.05	0.00	0.00	1.92	12.46	382.06	▇▁▁▁▁
yield_bu	52	0.68	59.21	11.57	30.44	49.23	62.93	68.64	74.80	▁▃▃▃▇
phyto_14	85	0.47	1.93	2.81	0.00	0.22	1.00	3.00	12.00	▇▁▁▁▁
dicambadrift_14	120	0.25	2.15	4.45	0.00	0.00	0.00	3.00	25.00	▇▁▁▁▁
canopyclosure_14	120	0.25	63.38	12.69	20.00	60.00	70.00	70.00	77.00	▁▁▁▅▇

Data visualization

Waterhemp control at 14 DAT

ggplot(data, aes(x=reorder(other,waterhempcontrol_14), y=waterhempcontrol_14,
       fill=trade, color=trade)) +
  geom_boxplot(color="black") +
  geom_jitter(alpha=0.2) +
  facet_grid(year ~ location) +
  coord_flip() +
  labs(x="", y="Waterhemp control (%)") +
  theme_minimal() +
  theme(legend.position = "none")

Waterhemp control at 28 DAT

ggplot(data, aes(x=reorder(other,waterhempcontrol_28), y=waterhempcontrol_28,
       fill=trade, color=trade)) +
  geom_boxplot(color="black") +
  geom_jitter(alpha=0.2) +
  facet_grid(year ~ location) +
  coord_flip() +
  labs(x="", y="Waterhemp control (%)") +
  theme_minimal() +
  theme(legend.position = "none")

Waterhemp control at harvest

ggplot(data, aes(x=reorder(other,waterhempcontrol_harvest), y=waterhempcontrol_harvest,
       fill=trade, color=trade)) +
  geom_boxplot(color="black") +
  geom_jitter(alpha=0.2) +
  facet_grid(year ~ location) +
  coord_flip() +
  labs(x="", y="Waterhemp control (%)") +
  theme_minimal() +
  theme(legend.position = "none")

Weed biomass

ggplot(data, aes(x=reorder(other,biomass_gm2), y=biomass_gm2,
       fill=trade, color=trade)) +
  geom_boxplot(color="black") +
  geom_jitter(alpha=0.2) +
  facet_grid(year ~ location) +
  coord_flip() +
  labs(x="", y="Weed biomass (g m2)") +
  theme_minimal() +
  theme(legend.position = "none")

Yield (bu/acre)

ggplot(data, aes(x=reorder(other,yield_bu), y=yield_bu,
       fill=trade, color=trade)) +
  geom_boxplot(color="black") +
  geom_jitter(alpha=0.2) +
  facet_grid(year ~ location) +
  coord_flip() +
  labs(x="", y="Waterhemp control (%)") +
  theme_minimal() +
  theme(legend.position = "none")

Data wrangling

new_dt <- 
  data %>% 
  rename (herbicide = other) %>% 
  mutate(
    wt_14 = waterhempcontrol_14/100,
    wt_28 = waterhempcontrol_28/100,
    yield_kg = yield_bu * 67.5) %>% 
  mutate(
    wt_14 = 
      case_when(
        waterhempcontrol_14 == 100   ~  0.99,
        TRUE            ~ wt_14),
    wt_28 = case_when(
        waterhempcontrol_28 == 100   ~  0.99,
        waterhempcontrol_28 == 0   ~  0.01,
        TRUE            ~ wt_28)) %>% 
  filter(
    herbicide != "PRE only"
  )

Data analysis

Waterhemp control at 14 DAT

Model

Model using herbicide trts, year, and location as fixed effects. Only rep as random effects

model_14 <- glmmTMB(wt_14 ~ herbicide * location * year + (1|rep), beta_family(link = "logit"), data=new_dt)

Anova

glmmTMB:::Anova.glmmTMB(model_14)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: wt_14
                            Chisq Df Pr(>Chisq)    
herbicide                 15.4841  8    0.05039 .  
location                   1.9251  1    0.16530    
year                    9714.1367  1    < 2e-16 ***
herbicide:location         6.4935  8    0.59213    
herbicide:year             1.4430  8    0.99361    
location:year              0.3624  1    0.54716    
herbicide:location:year    0.6957  8    0.99954    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova shows no interaction of year, herbicide and location. There is only herbicide and year effects

New model

Model using herbicide trt as fixed, and rep, year, and location as random effects

model_14_2 <- glmmTMB(wt_14 ~ herbicide + (1|rep) + (1|location/year), beta_family(link = "logit"), data=new_dt)

New Anova

glmmTMB:::Anova.glmmTMB(model_14_2)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: wt_14
          Chisq Df Pr(>Chisq)
herbicide 11.21  8     0.1901

No effects of herbicide treatments. However, I am looking into herbicide effects.

Least square means

lsmeans_14<- emmeans(model_14_2, ~ herbicide, cont="pairwise", adjust="none", type="response", alpha=0.05)

Compact letter display

plot(lsmeans_14, ~ herbicide, comparisons=TRUE, type="response", alpha=0.05, adjust="none")

CLD(lsmeans_14$emmeans, alpha=0.05, Letters=letters, adjust="none", reversed = TRUE) %>%
  kbl() %>%
  kable_classic_2(full_width = F)

	herbicide	response	SE	df	lower.CL	upper.CL	.group
3	PRE fb glufosinate + fomesafen + acetochlor	0.9663473	0.0092718	131	0.9423191	0.9805724	a
9	PRE fb glufosinate + fomesafen	0.9644099	0.0096027	131	0.9396867	0.9792228	a
1	PRE fb glufosinate + fomesafen + S-metolachlor	0.9633468	0.0100886	131	0.9372445	0.9788374	a
2	PRE fb glufosinate + pyroxasulfone	0.9513317	0.0122489	131	0.9205233	0.9705792	ab
7	PRE fb glufosinate + imazethapyr	0.9491984	0.0129393	131	0.9165874	0.9694842	ab
8	PRE fb glufosinate + acetochlor	0.9470465	0.0131010	131	0.9142870	0.9677273	ab
6	PRE fb glufosinate	0.9426195	0.0139116	131	0.9080580	0.9646943	ab
4	PRE fb glufosinate + S-metolachlor	0.9408544	0.0142213	131	0.9056144	0.9634678	ab
5	PRE fb glufosinate + dimethenamid-P	0.9296513	0.0162142	131	0.8900147	0.9557144	b

Waterhemp control at 28 DAT

Model

Model using herbicide trts, year, and location as fixed effects. Only rep as random effects

model_28 <- glmmTMB(wt_28 ~ herbicide * location * year + (1|rep), beta_family(link = "logit"), data=new_dt)

Anova

glmmTMB:::Anova.glmmTMB(model_28)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: wt_28
                           Chisq Df Pr(>Chisq)    
herbicide                23.1383 10   0.010247 *  
location                 14.2304  3   0.002608 ** 
year                    194.8150  1  < 2.2e-16 ***
herbicide:location       13.3607  8   0.100027    
herbicide:year            2.3217  8   0.969526    
location:year             0.0097  1   0.921564    
herbicide:location:year   0.0139  8   1.000000    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova shows no interaction of year, herbicide and location. There is only herbicide and year effects

New model

Model using herbicide trt as fixed, and rep, year, and location as random effects

model_28_2 <- glmmTMB(wt_28 ~ herbicide + (1|rep) + (1|location/year), beta_family(link = "logit"), data=new_dt)

New Anova

glmmTMB:::Anova.glmmTMB(model_28_2)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: wt_28
          Chisq Df Pr(>Chisq)    
herbicide 33.38  8  5.259e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No effects of herbicide treatments

Least square means

lsmeans_28 <- emmeans(model_28_2, ~ herbicide, cont="pairwise", adjust="none", type="response", alpha=0.05)

Plot

plot(lsmeans_28, ~ herbicide, comparisons=TRUE, type="response", alpha=0.05, adjust="none")

Compact letter display

CLD(lsmeans_28$emmeans, alpha=0.05, Letters=letters, adjust="none", reversed = TRUE) %>%
  kbl() %>%
  kable_classic_2(full_width = F)

	herbicide	response	SE	df	lower.CL	upper.CL	.group
1	PRE fb glufosinate + fomesafen	0.9335527	0.0219066	131	0.8747880	0.9658156	a
3	PRE fb glufosinate + fomesafen + S-metolachlor	0.9316031	0.0228235	131	0.8702296	0.9651143	a
7	PRE fb glufosinate + fomesafen + acetochlor	0.9297294	0.0230599	131	0.8681042	0.9637635	a
9	PRE fb glufosinate + acetochlor	0.9012060	0.0304182	131	0.8227128	0.9471779	ab
8	PRE fb glufosinate + pyroxasulfone	0.8942952	0.0319956	131	0.8124210	0.9429429	ab
2	PRE fb glufosinate + S-metolachlor	0.8814474	0.0351589	131	0.7925916	0.9353419	ab
5	PRE fb glufosinate + imazethapyr	0.8442439	0.0431294	131	0.7391023	0.9120553	bc
6	PRE fb glufosinate + dimethenamid-P	0.8436929	0.0424061	131	0.7407439	0.9106904	bc
4	PRE fb glufosinate	0.7817276	0.0528757	131	0.6598743	0.8686175	c

Yield

Yield raw data distribution

library(ggridges)
new_dt %>% 
  ggplot(
    aes(y=herbicide, x=yield_kg, fill=herbicide)) +
      geom_density_ridges(scale=2, show.legend = FALSE)

Data manipulation

yield_dt <- new_dt %>% 
  mutate(
    siteyr = case_when(
      location == "Lancaster" & year == 2020 ~ "Lan20",
      location == "Lancaster" & year == 2019 ~ "Lan19",
      location == "Brooklyn" & year == 2020 ~ "Bro20",
      location == "Brooklyn" & year == 2019 ~ "Bro19",
      TRUE                     ~  "siteyr"
    )
  )

Model

fit <- lmer(yield_kg ~ herbicide * siteyr + (1|rep), data=yield_dt)

Anova

anova(fit)

Type III Analysis of Variance Table with Satterthwaite's method
                   Sum Sq  Mean Sq NumDF  DenDF  F value Pr(>F)    
herbicide          428797    53600     8 67.094   0.5795 0.7912    
siteyr           35541129 17770564     2 67.950 192.1351 <2e-16 ***
herbicide:siteyr  1318370    82398    16 67.087   0.8909 0.5820    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Least square means

lsmeans_fit <- emmeans(fit, ~ siteyr, cont="pairwise", adjust="none", alpha=0.05)

Plot

plot(lsmeans_fit, ~ siteyr, comparisons=TRUE, type="response", alpha=0.05, adjust="none")

Compact letter display

CLD(lsmeans_fit$emmeans, alpha=0.05, Letters=letters, adjust="none", reversed = TRUE) %>%
  kbl() %>%
  kable_classic_2(full_width = F)

	siteyr	emmean	SE	df	lower.CL	upper.CL	.group
2	Lan20	4672.755	84.84491	4.867462	4452.856	4892.654	a
1	Bro19	4348.885	94.32107	6.946555	4125.503	4572.267	b
3	Bro20	3311.142	84.84491	4.867462	3091.242	3531.041	c

Soybean Glufosinate Study

Soybean Glufosinate Study

Load packages

Load dataset

Glimpse dataset

Skim dataset

Data visualization

Waterhemp control at 14 DAT

Waterhemp control at 28 DAT

Waterhemp control at harvest

Weed biomass

Yield (bu/acre)

Data wrangling

Data analysis

Waterhemp control at 14 DAT

Model

Anova

New model

New Anova

Least square means

Compact letter display

Waterhemp control at 28 DAT

Model

Anova

New model

New Anova

Least square means

Plot

Compact letter display

Yield

Yield raw data distribution

Data manipulation

Model

Anova

Least square means

Plot

Compact letter display